Load Packages

library(tictoc)
library(tidyverse)
library(cowplot)
library(ggExtra)
library(ggthemes)
library(ggrepel)
library(splatter)
library(scater)
library(scran)
library(kableExtra)
library(grid)
library(gridExtra)
library(r.jive)
library(Seurat)
library(harmony)

library(RSpectra)
library(scRNAseq)
library(lme4)

library(kBET)    # kBET
library(cluster) # ASW
library(lisi)    # LISI

source("PVCA.R")
source("JIVE/jive_speedup.R", chdir = T)

theme_set(theme_few())

dataset_name <- "dataset2"

Bacher T-Cell Data

# scRNAseq package
BacherTCellData <- BacherTCellData()

Cell Frequencies by Batch and Cell Cluster

batch_freq <- as.data.frame(table(Batch = BacherTCellData$batch))
batch_cluster_freq <-
  as.data.frame(table(
    Batch = BacherTCellData$batch,
    Cluster = BacherTCellData$new_cluster_names
  )) %>%
  pivot_wider(id_cols = c("Batch"), names_from = "Cluster", values_from = "Freq")

batch_cluster_df <- left_join(batch_freq, batch_cluster_freq, by = "Batch")
batch_cluster_df %>%
  kbl(
    col.names = c(
      "Batch Number",
      "Total Batch Count",
      "Central memory",
      "Cycling",
      "Cytotoxic / Th1",
      "Tfh-like",
      "Transitional memory",
      "Type-1 IFN signature"
    )
  ) %>%
  kable_styling() %>%
  add_header_above(c("Batch" = 2, "Cluster" = 6)) %>%
  column_spec(2, border_right = T)

Batch		Cluster
Batch Number	Total Batch Count	Central memory	Cycling	Cytotoxic / Th1	Tfh-like	Transitional memory	Type-1 IFN signature
1	2547	407	11	272	1399	445	13
2	11585	2027	261	1925	4438	2836	98
3	13344	1517	56	1894	6874	2815	188
4	12904	3324	24	1534	4608	3198	216
5	3059	955	13	331	921	736	103
10	587	47	1	326	130	80	3
11	7740	1196	11	1246	3479	1757	51
12	9961	1620	31	1373	4625	2255	57
13	4623	624	16	412	2243	1266	62
14	754	230	1	76	250	190	7
15	1113	239	4	143	425	289	13
16	376	87	0	61	138	87	3
17	296	65	0	41	107	82	1
18	17915	3275	125	2092	7954	4141	328
19	17613	2778	52	2536	8141	3988	118

We will test the integration methods on batches 14 and 15.

Preprocess Data

# Create Seurat object list so all data is preprocessed in the same manner
bacher_metadata <- as.data.frame(colData(BacherTCellData))
bacher_seurat <- CreateSeuratObject(counts(BacherTCellData), meta.data = bacher_metadata)

#################
# Identify HVGs #
#################

# Use Seurat feature selection method for selecting highly variable genes

# split the dataset into a list of seurat objects (one for each batch)
bacher_seurat_batch <- SplitObject(bacher_seurat, split.by = "batch")

# normalize and identify variable features for each dataset independently
bacher_seurat_batch <- lapply(X = bacher_seurat_batch, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = bacher_seurat_batch)

#######################################
# Create SCE object for visualization #
#######################################

# Select only HVGs
BacherTCellData <- BacherTCellData[features, ]
# Subset to batches of interest
BacherTCellData <- BacherTCellData[, colData(BacherTCellData)$batch %in% c(14, 15)]

# Final SCE object
data_SCE <- BacherTCellData
colData(data_SCE)$Batch <- factor(colData(data_SCE)$batch)
colData(data_SCE)$Cluster <- factor(colData(data_SCE)$new_cluster_names)

Visualize Raw Counts

data_SCE <- logNormCounts(data_SCE)

# Dimension reduction visualization
set.seed(1)
data_SCE <- runPCA(data_SCE, ncomponents = 30, ntop = 2000)
data_SCE <- runTSNE(data_SCE)
data_SCE <- runUMAP(data_SCE)

orig_title <- labs(title = "Original Data")

orig_p1 <- plot_grid(
  plotPCA(data_SCE, color_by = "Batch")   + scale_color_fivethirtyeight(name = "Batch") +
    orig_title,
  plotPCA(data_SCE, color_by = "Cluster") + scale_color_gdocs(name = "Cluster"),
  nrow=2, align = "v"
  ) + labs(title = "Test")

orig_p2 <- plot_grid(
  plotTSNE(data_SCE, color_by = "Batch")   + scale_color_fivethirtyeight(name = "Batch") +
    orig_title,
  plotTSNE(data_SCE, color_by = "Cluster") + scale_color_gdocs(name = "Cluster"),
  nrow=2, align = "v"
  )

orig_p3 <- plot_grid(
  plotUMAP(data_SCE, color_by = "Batch")   + scale_color_fivethirtyeight(name = "Batch") +
    orig_title,
  plotUMAP(data_SCE, color_by = "Cluster") + scale_color_gdocs(name = "Cluster"),
  nrow=2, align = "v"
  )

orig_p1 # ggsave("orig_pca.png", orig_p1)

orig_p2 # ggsave("orig_tsne.png", orig_p2)

orig_p3 # ggsave("orig_umap.png", orig_p3)

PVCA

batches <- data_SCE$Batch
clusters <- data_SCE$Cluster
rawcounts <- as.matrix(logcounts(data_SCE))

# Perform PVCA for 10 random samples of size 1000 (more computationally efficient)
pvca_res <- matrix(nrow = 10, ncol = 3)
meta <- cbind(Batch = batches, Cluster = clusters)

set.seed(1)
for (i in 1:10){
  sample <- sample(1:ncol(rawcounts), 1000, replace = FALSE)
  pvca_res[i, ] <- PVCA(rawcounts[, sample], meta[sample, ], threshold = 0.6, inter = FALSE)
}

# Average effect size across samples
pvca_means <- colMeans(pvca_res)
names(pvca_means) <- c(colnames(meta), "Residual")

# Plot PVCA
pvca_plot <- PlotPVCA(pvca_means, "PVCA of Raw Count Data")
pvca_plot

Perform Integration Methods

JIVE

method_name <- "JIVE"

Run Integration Method

n_batches <- length(unique(batches))
unq_batches <- sort(unique(batches))

jive_data <- NULL
all_clusters <- NULL

for (i in 1:n_batches) {
  jive_data[[as.character(unq_batches[i])]] <- t(rawcounts[, batches == unq_batches[i]])
  all_clusters[[as.character(unq_batches[i])]] <- clusters[batches == unq_batches[i]]
}

CORES <- parallel::detectCores()

tic("JIVE v2 runtime (given)")
JIVE_results <-
  jive_v2(
    jive_data,
    rankJ = 5,
    rankA = rep(10, length(jive_data)),
    method = "given",
    maxiter = 10000,
    CORES = CORES
  )
# JIVE v2 runtime (given): 277.38 sec elapsed
# JIVE algorithm converged after  245  iterations
toc()

saveRDS(JIVE_results, file = "JIVE/data/JIVE_v2_dataset2_b14_b15.rds")

Summary of estimated ranks and variance explained by joint/individual/residual for each approach:

JIVE_results <- readRDS(file = "JIVE/data/JIVE_v2_dataset2_b14_b15.rds")
summary(JIVE_results)

## $Method
## [1] "given"
## 
## $Ranks
##      Source  Rank
## [1,] "Joint" "5" 
## [2,] "14"    "10"
## [3,] "15"    "10"
## 
## $Variance
##               14    15
## Joint      0.541 0.549
## Individual 0.080 0.040
## Residual   0.380 0.410

plot(JIVE_results)

joint <- do.call(rbind, JIVE_results$joint)
jive_final <- joint

Visualization of Integration

joint_svd <- svds(joint, k = JIVE_results$rankJ)

var_explained <- data.frame(
  sv = joint_svd$d,
  cum_sv = cumsum(joint_svd$d),
  tot_sv = sum(joint_svd$d)
  ) %>%
  mutate(
    pct_sv = sv / tot_sv,
    cumpct_sv = cum_sv / tot_sv,
    PC = paste0("PC", row_number()),
    PC = factor(PC, levels = paste0("PC", row_number()))
  )

var_explained %>%
  ggplot(aes(x = PC)) +
  geom_col(aes(y = cumpct_sv, fill = PC), color = "black") +
  geom_col(aes(y = cumpct_sv - pct_sv), fill = "gray", color = "black") +
  geom_label(aes(y = cumpct_sv - (pct_sv / 2), label = paste0(round(pct_sv * 100, 2), "%"))) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    x = "", y = "",
    title = "Cumulative Percentage of Total Variance Explained",
    subtitle = "Integrated (Joint) Data"
  )

jive_integrated <- SingleCellExperiment(list(logcounts = t(jive_final)))
colData(jive_integrated)$Batch <- batches
colData(jive_integrated)$Cluster <- clusters

set.seed(1)
# Original PCA based on joint matrices for each batch
jive_integrated <- runPCA(jive_integrated, ncomponents = 30, ntop = 2000)

# t-SNE based on original PC
jive_integrated <- runTSNE(jive_integrated)

# UMAP based on original PC
jive_integrated <- runUMAP(jive_integrated)

plot_title <- labs(
    title = "Integrated Data",
    subtitle = paste0(
      "Joint Rank: ",
      JIVE_results$rankJ,
      ", Individual Ranks: ",
      paste(JIVE_results$rankA, collapse = ", ")
    )
  )

p1_options <- scale_color_fivethirtyeight(name = "Batch")
p2_options <- scale_color_gdocs(name = "Cluster")

# PCA
jive_pca <- plot_grid(
  plotPCA(jive_integrated, color_by = "Batch")   + plot_title + p1_options,
  plotPCA(jive_integrated, color_by = "Cluster") + p2_options,
  nrow = 2, align = "v"
  )

# jive_pca_2_3 <- plot_grid(
#   plotPCA(jive_integrated, color_by = "Batch",   ncomponents = 2:3) + plot_title + p1_options,
#   plotPCA(jive_integrated, color_by = "Cluster", ncomponents = 2:3) + p2_options,
#   nrow = 2, align = "v"
#   )
# 
# jive_pca_3_4 <- plot_grid(
#   plotPCA(jive_integrated, color_by = "Batch",   ncomponents = 3:4) + plot_title + p1_options,
#   plotPCA(jive_integrated, color_by = "Cluster", ncomponents = 3:4) + p2_options,
#   nrow = 2, align = "v"
#   )
# 
# jive_pca_4_5 <- plot_grid(
#   plotPCA(jive_integrated, color_by = "Batch",   ncomponents = 4:5) + plot_title + p1_options,
#   plotPCA(jive_integrated, color_by = "Cluster", ncomponents = 4:5) + p2_options,
#   nrow = 2, align = "v"
#   )

# t-SNE
jive_tsne <- plot_grid(
  plotTSNE(jive_integrated, color_by = "Batch")   + p1_options,
  plotTSNE(jive_integrated, color_by = "Cluster") + p2_options,
  nrow = 2, align = "v"
  )

# UMAP
jive_umap <- plot_grid(
  plotUMAP(jive_integrated, color_by = "Batch")   + p1_options,
  plotUMAP(jive_integrated, color_by = "Cluster") + p2_options,
  nrow = 2, align = "v"
  )

jive_pca

jive_tsne

jive_umap

kBET

Sub-sample data at varying percentages
- 5%, 10%, 15%, 20%, 25%
- Stratified sampling so that each batch has the same number of cells selected
  - E.g., assume 2 batches, sampling 10% of 1000 cells = 100 cells: sample 50 cells from each batch
Calculate kBET rejection rates at each percentage
- Lower RR indicate well-mixed batches
- Higher RR indicate poorly-mixed batches

# data: a matrix (rows: samples, columns: features (genes))
data <- svds(jive_final, k = 30)$u

# batch: vector or factor with batch label of each cell 
batch <- batches

sample_size <- seq(0.05, 0.25, by = 0.05)
rejection_rate <- list()

set.seed(1)
for (i in 1:length(sample_size)) {
  k <- floor(nrow(data) * sample_size[i])
  
  cat("\nSample Size: ", sample_size[i], ", k = ", k, "\n", sep = "")
  cat("------------------------------\n")
 
  batch.estimate <- kBET(data, batch, k0 = k, plot = F, do.pca = F, verbose = T)
  
  cat("\n------------------------------\n")
  
  if (class(batch.estimate) == "list") {
    rejection_rate[[i]] <- batch.estimate$summary %>%
      rownames_to_column(var = "statistic") %>%
      mutate(sample_size = sample_size[i]) %>%
      select(sample_size, statistic, kBET_rr = kBET.observed)
  }
}

## 
## Sample Size: 0.05, k = 93
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.23    0.00    0.23 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.1, k = 186
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.33    0.00    0.33 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.15, k = 280
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.41    0.00    0.41 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.2, k = 373
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.51    0.00    0.52 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.25, k = 466
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.63    0.01    0.64 
## Number of kBET tests is set to 187.
## 
## ------------------------------

jive_kbet <- bind_rows(rejection_rate) %>%
  mutate(Method = method_name) %>%
  pivot_wider(
    id_cols = c("Method", "sample_size"),
    names_from = "statistic",
    values_from = "kBET_rr"
  ) %>%
  select(Method, sample_size, median_rr = `50%`)

jive_kbet %>% write_csv(file = paste0("output/", method_name, "_", dataset_name, "_kBET.csv"))

Average Silhouette Width

Original batch label and cell types/groups will be considered the clusters to compute \(ASW_{batch}\) and \(ASW_{group}\)

Compute average silhouette width (ASW) for the first 30 PCs of 80% sub-sample of joint data
- ASW ranges between -1 and 1
- Values close to 1 indicate the cell is well-matched to its cluster (i.e., batch label)
- Values close to -1 indicate the cell is not well-matched to its cluster (i.e., batch label)
  - We hope to see lower values for \(ASW_{batch}\) since it is indicative of well-mixed batches
  - We hope to see higher values for \(ASW_{group}\) since it is indicative that distinct cell groups were preserved after batch-mixing
Repeat 20 times
Calculate the median ASW from the 20 runs to ensure stability of the measurement

asw_repeat <- 20
npcs <- 30

asw_batch <- list()
asw_cluster <- list()
subset_prop <- 0.8 # subsample to 80% of the data
subset_size_total <- floor(length(batch) * subset_prop)

for (i in 1:asw_repeat) {
  set.seed(i)
  subset_index <- 
    data.frame(index = 1:length(batch) , batch = batch) %>%
    slice_sample(n = subset_size_total, replace = F)
  
  subset_id <- subset_index %>%
    pull(index)
  
  data_pc <- svds(jive_final[subset_id, ], k = npcs)$u
  
  dissimilarity_matrix <- daisy(data_pc)
  
  asw_batch_i <- silhouette(as.integer(factor(batches[subset_id])), dissimilarity_matrix)
  asw_batch[[i]] <- as.data.frame(rbind(summary(asw_batch_i)[["si.summary"]]))
  
  asw_cluster_i <- silhouette(as.integer(factor(clusters[subset_id])), dissimilarity_matrix)
  asw_cluster[[i]] <- as.data.frame(rbind(summary(asw_cluster_i)[["si.summary"]]))
}

jive_asw <- bind_rows(
  bind_rows(asw_batch) %>% mutate(Label = "Batch"),
  bind_rows(asw_cluster) %>% mutate(Label = "Cluster")
  ) %>%
  mutate(Method = method_name) %>%
  select(Method, Label, Mean)

jive_asw %>% write_csv(file = paste0("output/", method_name, "_", dataset_name, "_ASW.csv"))

LISI Scores

# Performed for all cells
jive_lisi <-
  compute_lisi(
    jive_final,
    data.frame(Batch = batches, Cluster = clusters),
    c("Batch", "Cluster"),
    perplexity = 40
  ) %>%
  pivot_longer(
    cols = c("Batch", "Cluster"),
    names_to = "Label",
    values_to = "LISI"
  ) %>%
  mutate(Method = method_name)

jive_lisi %>%
  write_csv(file = paste0("output/", method_name, "_", dataset_name, "_LISI.csv"))

Seurat v3

method_name <- "Seurat"

batch_indices <- which(unique(bacher_seurat$batch) %in% c(14, 15)) # subset batches
bacher_seurat_batch <- bacher_seurat_batch[batch_indices]

Run Integration Method

Perform integration by finding a set of anchors from the two batches (via CCA)
Use the anchor list to integrate the two batches

anchors <- FindIntegrationAnchors(object.list = bacher_seurat_batch, anchor.features = features)

## Scaling features for provided objects

## Finding all pairwise anchors

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 3141 anchors

## Filtering anchors

##  Retained 1835 anchors

seurat_integrated <- IntegrateData(anchorset = anchors)

## Merging dataset 1 into 2

## Extracting anchors for merged samples

## Finding integration vectors

## Finding integration vector weights

## Integrating data

GetAssay(seurat_integrated, "integrated")

## Assay data with 2000 features for 1867 cells
## Top 10 variable features:
##  CCL4, IL3, CCL4L2, XCL1, XCL2, CSF2, IL10, CCL3, GZMB, IL4

Integrated/combined data is centered/scaled
First 30 PCs are computed/stored
PC1-PC30 are used to perform t-SNE and UMAP dimensionality reduction
Calculates k-nearest neighbors (k = 20 by default) from PC1-PC30
Identify clusters by using Shared Nearest Neighbors (SNN):
- “First calculate k-nearest neighbors and construct the SNN graph”
- “Then optimize the modularity function to determine clusters”

# specify that we will perform downstream analysis on the corrected data note that the
# original unmodified data still resides in the 'RNA' assay
DefaultAssay(seurat_integrated) <- "integrated"

# Run the standard workflow for visualization and clustering
set.seed(1)
seurat_integrated <- ScaleData(seurat_integrated, verbose = FALSE)
seurat_integrated <- RunPCA(seurat_integrated, verbose = FALSE, approx = FALSE)
seurat_integrated <- RunTSNE(seurat_integrated, reduction = "pca", dims = 1:30)
seurat_integrated <- RunUMAP(seurat_integrated, reduction = "pca", dims = 1:30)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 09:03:56 UMAP embedding parameters a = 0.9922 b = 1.112

## 09:03:56 Read 1867 rows and found 30 numeric columns

## 09:03:56 Using Annoy for neighbor search, n_neighbors = 30

## 09:03:56 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 09:03:56 Writing NN index file to temp file C:\Users\Joseph\AppData\Local\Temp\RtmpKKdqng\file27385ec14025
## 09:03:56 Searching Annoy index using 1 thread, search_k = 3000
## 09:03:57 Annoy recall = 100%
## 09:03:58 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:03:59 Initializing from normalized Laplacian + noise (using RSpectra)
## 09:03:59 Commencing optimization for 500 epochs, with 75500 positive edges
## 09:04:04 Optimization finished

seurat_integrated <- FindNeighbors(seurat_integrated, reduction = "pca", dims = 1:30)

## Computing nearest neighbor graph
## Computing SNN

seurat_integrated <- FindClusters(seurat_integrated, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1867
## Number of edges: 98822
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7864
## Number of communities: 7
## Elapsed time: 0 seconds

seurat_integrated <- AddMetaData(seurat_integrated, factor(seurat_integrated$batch), col.name = "Batch")
seurat_integrated <- AddMetaData(seurat_integrated, factor(seurat_integrated$new_cluster_names), col.name = "Cluster")

seurat_final <- t(seurat_integrated@assays$integrated@scale.data)

Visualization of Integration

plot_title <- labs(
    title = "Integrated Data",
    subtitle = method_name
  )

# PCA
seurat_pca <- plot_grid(
  PCAPlot(seurat_integrated, group.by = "Batch")   + plot_title  + p1_options,
  PCAPlot(seurat_integrated, group.by = "Cluster") + ggtitle("") + p2_options,
  nrow = 2, align = "vh"
  )

# t-SNE
seurat_tsne <- plot_grid(
  TSNEPlot(seurat_integrated, group.by = "Batch")   + plot_title  + p1_options,
  TSNEPlot(seurat_integrated, group.by = "Cluster") + ggtitle("") + p2_options,
  nrow = 2, align = "vh"
  )

# UMAP
seurat_umap <- plot_grid(
  UMAPPlot(seurat_integrated, group.by = "Batch")   + plot_title  + p1_options,
  UMAPPlot(seurat_integrated, group.by = "Cluster") + ggtitle("") + p2_options,
  nrow = 2, align = "vh"
  )

seurat_pca

seurat_tsne

seurat_umap

kBET

# data: a matrix (rows: samples, columns: features (genes))
data <- svds(t(seurat_integrated@assays$integrated@scale.data), k = 30)$u

# batch: vector or factor with batch label of each cell 
batch <- batches

sample_size <- seq(0.05, 0.25, by = 0.05)
rejection_rate <- list()

set.seed(1)
for (i in 1:length(sample_size)) {
  k <- floor(nrow(data) * sample_size[i])
  
  cat("\nSample Size: ", sample_size[i], ", k = ", k, "\n", sep = "")
  cat("------------------------------\n")
 
  batch.estimate <- kBET(data, batch, k0 = k, plot = F, do.pca = F, verbose = T)
  
  cat("\n------------------------------\n")
  
  if (class(batch.estimate) == "list") {
    rejection_rate[[i]] <- batch.estimate$summary %>%
      rownames_to_column(var = "statistic") %>%
      mutate(sample_size = sample_size[i]) %>%
      select(sample_size, statistic, kBET_rr = kBET.observed)
  }
}

## 
## Sample Size: 0.05, k = 93
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.29    0.00    0.29 
## Number of kBET tests is set to 187.
## There are 24 cells (1.285%) that do not appear in any neighbourhood.
## The expected frequencies for each category have been adapted.
## Cell indexes are saved to result list.
## 
## ------------------------------
## 
## Sample Size: 0.1, k = 186
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.39    0.00    0.40 
## Number of kBET tests is set to 187.
## There are 13 cells (0.696%) that do not appear in any neighbourhood.
## The expected frequencies for each category have been adapted.
## Cell indexes are saved to result list.
## 
## ------------------------------
## 
## Sample Size: 0.15, k = 280
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.48    0.00    0.49 
## Number of kBET tests is set to 187.
## There are 11 cells (0.589%) that do not appear in any neighbourhood.
## The expected frequencies for each category have been adapted.
## Cell indexes are saved to result list.
## 
## ------------------------------
## 
## Sample Size: 0.2, k = 373
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.61    0.00    0.61 
## Number of kBET tests is set to 187.
## There are 9 cells (0.482%) that do not appear in any neighbourhood.
## The expected frequencies for each category have been adapted.
## Cell indexes are saved to result list.
## 
## ------------------------------
## 
## Sample Size: 0.25, k = 466
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.75    0.00    0.75 
## Number of kBET tests is set to 187.
## There are 7 cells (0.375%) that do not appear in any neighbourhood.
## The expected frequencies for each category have been adapted.
## Cell indexes are saved to result list.
## 
## ------------------------------

seurat_kbet <- bind_rows(rejection_rate) %>%
  mutate(Method = method_name) %>%
  pivot_wider(
    id_cols = c("Method", "sample_size"),
    names_from = "statistic",
    values_from = "kBET_rr"
  ) %>%
  select(Method, sample_size, median_rr = `50%`)

seurat_kbet %>% write_csv(file = paste0("output/", method_name, "_", dataset_name, "_kBET.csv"))

Average Silhouette Width

asw_repeat <- 20
npcs <- 30

asw_batch <- list()
asw_cluster <- list()
subset_prop <- 0.8 # subsample to 80% of the data
subset_size_total <- floor(length(batch) * subset_prop)

for (i in 1:asw_repeat) {
  set.seed(i)
  subset_index <- 
    data.frame(index = 1:length(batch) , batch = batch) %>%
    slice_sample(n = subset_size_total, replace = F)
  
  subset_id <- subset_index %>%
    pull(index)
  
  seurat_tmp <- NULL
  seurat_tmp <-
    RunPCA(
      seurat_integrated[, subset_id],
      npcs = npcs,
      verbose = FALSE,
      reduction.name = "subset_pca",
      reduction.key = "sPC_",
      seed.use = NULL
    )
  data_pc <- seurat_tmp@reductions$subset_pca@cell.embeddings
  
  dissimilarity_matrix <- daisy(data_pc)
  
  asw_batch_i <- silhouette(as.integer(factor(batches[subset_id])), dissimilarity_matrix)
  asw_batch[[i]] <- as.data.frame(rbind(summary(asw_batch_i)[["si.summary"]]))
  
  asw_cluster_i <- silhouette(as.integer(factor(clusters[subset_id])), dissimilarity_matrix)
  asw_cluster[[i]] <- as.data.frame(rbind(summary(asw_cluster_i)[["si.summary"]]))
}

seurat_asw <- bind_rows(
  bind_rows(asw_batch) %>% mutate(Label = "Batch"),
  bind_rows(asw_cluster) %>% mutate(Label = "Cluster")
  ) %>%
  mutate(Method = method_name) %>%
  select(Method, Label, Mean)

seurat_asw %>% write_csv(file = paste0("output/", method_name, "_", dataset_name, "_ASW.csv"))

LISI Scores

# Performed for all cells
seurat_lisi <-
  compute_lisi(
    seurat_final,
    data.frame(Batch = batches, Cluster = clusters),
    c("Batch", "Cluster"),
    perplexity = 40
  ) %>%
  pivot_longer(
    cols = c("Batch", "Cluster"),
    names_to = "Label",
    values_to = "LISI"
  ) %>%
  mutate(Method = method_name)

seurat_lisi %>%
  write_csv(file = paste0("output/", method_name, "_", dataset_name, "_LISI.csv"))

Harmony

method_name <- "Harmony"

Run Integration Method

harmony_integrated <-
  runPCA(
    data_SCE,
    ncomponents = ncol(data_SCE),
    ntop = nrow(data_SCE),
    BSPARAM = BiocSingular::ExactParam()
  )
harmony_integrated <-
  RunHarmony(
    harmony_integrated,
    "Batch",
    lambda = 0.1
  )

## Harmony 1/10

## Harmony 2/10

## Harmony 3/10

## Harmony 4/10

## Harmony 5/10

## Harmony 6/10

## Harmony converged after 6 iterations

harmony_integrated <- runTSNE(harmony_integrated, dimred = "HARMONY")
harmony_integrated <- runUMAP(harmony_integrated, dimred = "HARMONY")

harmony_final <- reducedDim(harmony_integrated, "HARMONY")

Visualization of Integration

plot_title <- labs(title = "Integrated Data", subtitle = method_name)

# Harmony PCA
harmony_pca <- plot_grid(
  plotReducedDim(harmony_integrated, dimred = "HARMONY",
    color_by = "Batch",   ncomponents = 1:2 ) + plot_title + p1_options,
  plotReducedDim(harmony_integrated, dimred = "HARMONY",
    color_by = "Cluster", ncomponents = 1:2 ) + p2_options,
  nrow = 2, align = "vh"
  )

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

# t-SNE
harmony_tsne <- plot_grid(
  plotTSNE(harmony_integrated, color_by = "Batch")   + plot_title + p1_options,
  plotTSNE(harmony_integrated, color_by = "Cluster") + p2_options,
  nrow = 2, align = "vh"
  )

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

# UMAP
harmony_umap <- plot_grid(
  plotUMAP(harmony_integrated, color_by = "Batch")   + plot_title + p1_options,
  plotUMAP(harmony_integrated, color_by = "Cluster") + p2_options,
  nrow = 2, align = "vh"
  )

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

harmony_pca

harmony_tsne

harmony_umap

kBET

# data: a matrix (rows: samples, columns: features (genes))
data <- harmony_final[, 1:30]

# batch: vector or factor with batch label of each cell 
batch <- batches

sample_size <- seq(0.05, 0.25, by = 0.05)
rejection_rate <- list()

set.seed(1)
for (i in 1:length(sample_size)) {
  k <- floor(nrow(data) * sample_size[i])
  
  cat("\nSample Size: ", sample_size[i], ", k = ", k, "\n", sep = "")
  cat("------------------------------\n")
 
  batch.estimate <- kBET(data, batch, k0 = k, plot = F, do.pca = F, verbose = T)
  
  cat("\n------------------------------\n")
  
  if (class(batch.estimate) == "list") {
    rejection_rate[[i]] <- batch.estimate$summary %>%
      rownames_to_column(var = "statistic") %>%
      mutate(sample_size = sample_size[i]) %>%
      select(sample_size, statistic, kBET_rr = kBET.observed)
  }
}

## 
## Sample Size: 0.05, k = 93
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.25    0.00    0.25 
## Number of kBET tests is set to 187.
## No outsiders found.
## ------------------------------
## 
## Sample Size: 0.1, k = 186
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.32    0.00    0.31 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.15, k = 280
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.39    0.00    0.39 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.2, k = 373
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.48    0.02    0.50 
## Number of kBET tests is set to 187.
## 
## ------------------------------
## 
## Sample Size: 0.25, k = 466
## ------------------------------
## finding knns...done. Time:
##    user  system elapsed 
##    0.58    0.01    0.59 
## Number of kBET tests is set to 187.
## 
## ------------------------------

harmony_kbet <- bind_rows(rejection_rate) %>%
  mutate(Method = method_name) %>%
  pivot_wider(
    id_cols = c("Method", "sample_size"),
    names_from = "statistic",
    values_from = "kBET_rr"
  ) %>%
  select(Method, sample_size, median_rr = `50%`)

harmony_kbet %>% write_csv(file = paste0("output/", method_name, "_", dataset_name, "_kBET.csv"))

Average Silhouette Width

asw_repeat <- 20
npcs <- 30

asw_batch <- list()
asw_cluster <- list()
subset_prop <- 0.8 # subsample to 80% of the data
subset_size_total <- floor(length(batch) * subset_prop)

for (i in 1:asw_repeat) {
  set.seed(i)
  subset_index <- 
    data.frame(index = 1:length(batch) , batch = batch) %>%
    slice_sample(n = subset_size_total, replace = F)
  
  subset_id <- subset_index %>%
    pull(index)
  
  data_pc <- svds(harmony_final[subset_id, ], k = npcs)$u
  
  dissimilarity_matrix <- daisy(data_pc)
  
  asw_batch_i <- silhouette(as.integer(factor(batches[subset_id])), dissimilarity_matrix)
  asw_batch[[i]] <- as.data.frame(rbind(summary(asw_batch_i)[["si.summary"]]))
  
  asw_cluster_i <- silhouette(as.integer(factor(clusters[subset_id])), dissimilarity_matrix)
  asw_cluster[[i]] <- as.data.frame(rbind(summary(asw_cluster_i)[["si.summary"]]))
}

harmony_asw <- bind_rows(
  bind_rows(asw_batch) %>% mutate(Label = "Batch"),
  bind_rows(asw_cluster) %>% mutate(Label = "Cluster")
  ) %>%
  mutate(Method = method_name) %>%
  select(Method, Label, Mean)

harmony_asw %>% write_csv(file = paste0("output/", method_name, "_", dataset_name, "_ASW.csv"))

LISI Scores

# Performed for all cells
harmony_lisi <-
  compute_lisi(
    harmony_final,
    data.frame(Batch = batches, Cluster = clusters),
    c("Batch", "Cluster"),
    perplexity = 40
  ) %>%
  pivot_longer(
    cols = c("Batch", "Cluster"),
    names_to = "Label",
    values_to = "LISI"
  ) %>%
  mutate(Method = method_name)

harmony_lisi %>%
  write_csv(file = paste0("output/", method_name, "_", dataset_name, "_LISI.csv"))

Integration Method Comparison

method_lvl <- c("JIVE", "Seurat", "Harmony", "Raw")

t-SNE Plots

# JIVE
jive_tsne_df <- reducedDim(jive_integrated, "TSNE") %>%
  as.data.frame() %>%
  mutate(Method = "JIVE") %>%
  mutate(Batch = batches, Cluster = clusters)

jive_tsne_batch <- jive_tsne_df %>%
  ggplot(aes(x = V1, y = V2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), legend.position = "right") +
  labs(title = "JIVE") + p1_options

jive_tsne_cluster <- jive_tsne_df %>%
  ggplot(aes(x = V1, y = V2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(legend.position = "right") + p2_options

jive_tsne_plot <- plot_grid(
  jive_tsne_batch + theme(legend.position = "none"),
  jive_tsne_cluster + theme(legend.position = "none"),
  nrow = 2, align = "hv"
)

# Seurat
seurat_tsne_df <- Embeddings(seurat_integrated, "tsne") %>%
  as.data.frame() %>%
  mutate(Method = "Seurat") %>%
  mutate(Batch = batches, Cluster = clusters)

seurat_tsne_batch <- seurat_tsne_df %>%
  ggplot(aes(x = tSNE_1, y = tSNE_2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = "Seurat") + p1_options

seurat_tsne_cluster <- seurat_tsne_df %>%
  ggplot(aes(x = tSNE_1, y = tSNE_2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  p2_options

seurat_tsne_plot <- plot_grid(
  seurat_tsne_batch,
  seurat_tsne_cluster,
  nrow = 2, align = "hv"
)

# Harmony
harmony_tsne_df <- reducedDim(harmony_integrated, "TSNE") %>%
  as.data.frame() %>%
  mutate(Method = "Harmony") %>%
  mutate(Batch = batches, Cluster = clusters)

harmony_tsne_batch <- harmony_tsne_df %>%
  ggplot(aes(x = V1, y = V2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = "Harmony") + p1_options

harmony_tsne_cluster <- harmony_tsne_df %>%
  ggplot(aes(x = V1, y = V2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  p2_options

harmony_tsne_plot <- plot_grid(
  harmony_tsne_batch,
  harmony_tsne_cluster,
  nrow = 2, align = "hv"
)

# Raw
raw_tsne_df <- reducedDim(data_SCE, "TSNE") %>%
  as.data.frame() %>%
  mutate(Method = "Raw") %>%
  mutate(Batch = batches, Cluster = clusters)

raw_tsne_batch <- raw_tsne_df %>%
  ggplot(aes(x = V1, y = V2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = "Raw") + p1_options

raw_tsne_cluster <- raw_tsne_df %>%
  ggplot(aes(x = V1, y = V2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  p2_options

raw_tsne_plot <- plot_grid(
  raw_tsne_batch,
  raw_tsne_cluster,
  nrow = 2, align = "hv"
)

# All
legend_batch <- get_legend(
  # create some space to the left of the legend
  jive_tsne_batch + theme(legend.box.margin = margin(0, 0, 0, 12))
)

legend_cluster <- get_legend(
  # create some space to the left of the legend
  jive_tsne_cluster + theme(legend.box.margin = margin(0, 0, 0, 12))
)

tsne_prow <- plot_grid(
  jive_tsne_plot,
  seurat_tsne_plot,
  harmony_tsne_plot,
  raw_tsne_plot,
  nrow = 1,
  hjust = -1
  )

tsne_legend_both <- plot_grid(legend_batch, legend_cluster, nrow = 2, align = "v")

all_tsne_plot <- plot_grid(tsne_prow, tsne_legend_both, rel_widths = c(3, 1))

grid.arrange(arrangeGrob(
  all_tsne_plot,
  left = textGrob("t-SNE 2", rot = 90),
  bottom = textGrob("t-SNE 1")
))

UMAP Plots

# JIVE
jive_umap_df <- reducedDim(jive_integrated, "UMAP") %>%
  as.data.frame() %>%
  mutate(Method = "JIVE") %>%
  mutate(Batch = batches, Cluster = clusters)

jive_umap_batch <- jive_umap_df %>%
  ggplot(aes(x = V1, y = V2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), legend.position = "right") +
  labs(title = "JIVE") + p1_options

jive_umap_cluster <- jive_umap_df %>%
  ggplot(aes(x = V1, y = V2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(legend.position = "right") + p2_options

jive_umap_plot <- plot_grid(
  jive_umap_batch + theme(legend.position = "none"),
  jive_umap_cluster + theme(legend.position = "none"),
  nrow = 2, align = "hv"
)

# Seurat
seurat_umap_df <- Embeddings(seurat_integrated, "umap") %>%
  as.data.frame() %>%
  mutate(Method = "Seurat") %>%
  mutate(Batch = batches, Cluster = clusters)

seurat_umap_batch <- seurat_umap_df %>%
  ggplot(aes(x = UMAP_1, y = UMAP_2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = "Seurat") + p1_options

seurat_umap_cluster <- seurat_umap_df %>%
  ggplot(aes(x = UMAP_1, y = UMAP_2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  p2_options

seurat_umap_plot <- plot_grid(
  seurat_umap_batch,
  seurat_umap_cluster,
  nrow = 2, align = "hv"
)

# Harmony
harmony_umap_df <- reducedDim(harmony_integrated, "UMAP") %>%
  as.data.frame() %>%
  mutate(Method = "Harmony") %>%
  mutate(Batch = batches, Cluster = clusters)

harmony_umap_batch <- harmony_umap_df %>%
  ggplot(aes(x = V1, y = V2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = "Harmony") + p1_options

harmony_umap_cluster <- harmony_umap_df %>%
  ggplot(aes(x = V1, y = V2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  p2_options

harmony_umap_plot <- plot_grid(
  harmony_umap_batch,
  harmony_umap_cluster,
  nrow = 2, align = "hv"
)

# Raw
raw_umap_df <- reducedDim(data_SCE, "UMAP") %>%
  as.data.frame() %>%
  mutate(Method = "Raw") %>%
  mutate(Batch = batches, Cluster = clusters)

raw_umap_batch <- raw_umap_df %>%
  ggplot(aes(x = V1, y = V2, color = Batch)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = "Raw") + p1_options

raw_umap_cluster <- raw_umap_df %>%
  ggplot(aes(x = V1, y = V2, color = Cluster)) +
  geom_point(alpha = 0.5) +
  theme_nothing() +
  p2_options

raw_umap_plot <- plot_grid(
  raw_umap_batch,
  raw_umap_cluster,
  nrow = 2, align = "hv"
)

# All
legend_batch <- get_legend(
  # create some space to the left of the legend
  jive_umap_batch + theme(legend.box.margin = margin(0, 0, 0, 12))
)

legend_cluster <- get_legend(
  # create some space to the left of the legend
  jive_umap_cluster + theme(legend.box.margin = margin(0, 0, 0, 12))
)

umap_prow <- plot_grid(
  jive_umap_plot,
  seurat_umap_plot,
  harmony_umap_plot,
  raw_umap_plot,
  nrow = 1,
  hjust = -1
  )

umap_legend_both <- plot_grid(legend_batch, legend_cluster, nrow = 2, align = "v")

all_umap_plot <- plot_grid(umap_prow, umap_legend_both, rel_widths = c(3, 1))

grid.arrange(arrangeGrob(
  all_umap_plot,
  left = textGrob("UMAP 2", rot = 90),
  bottom = textGrob("UMAP 1")
))

kBET

all_kbet <- bind_rows(jive_kbet, seurat_kbet, harmony_kbet, raw_kbet) %>%
  mutate(Method = factor(Method, levels = method_lvl))

kbet_plot <- all_kbet %>%
  ggplot(aes(
    x = sample_size,
    y = I(1 - median_rr),
    group = Method,
    color = Method
  )) +  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_x_continuous(labels = scales::percent) +
  expand_limits(y = c(0, 1)) +
  labs(
    x = "% sample size",
    y = "kBET (acceptance rate)",
    color = "Integration Method"
    )

kbet_plot + labs(title = "kBET Acceptance Rates")

Average Silhouette Width (ASW)

all_asw <- bind_rows(jive_asw, seurat_asw, harmony_asw, raw_asw) %>%
  mutate(Method = factor(Method, levels = method_lvl))

asw_plot <- all_asw %>%
  group_by(Label) %>%
  mutate(Mean = scales::rescale(Mean)) %>%
  group_by(Method, Label) %>%
  summarize(median_ASW = median(Mean), .groups = "keep") %>%
  pivot_wider(
    id_cols = "Method",
    names_from = "Label",
    values_from = "median_ASW"
  ) %>% 
  ggplot(aes(x = Cluster, y = I(1 - Batch), color = Method)) +
  geom_point(size = 2) +
  geom_text_repel(aes(label = Method), show.legend = FALSE) +
  labs(x = "ASW cell type", y = "1 - ASW batch")

asw_plot + labs(title = "Average Silhouette Width Scores")

LISI Scores

all_lisi <- bind_rows(jive_lisi, seurat_lisi, harmony_lisi, raw_lisi) %>%
  mutate(Method = factor(Method, levels = method_lvl))

lisi_plot <- all_lisi %>%
  group_by(Label) %>%
  mutate(LISI = scales::rescale(LISI)) %>%
  group_by(Method, Label) %>%
  summarize(median_LISI = median(LISI), .groups = "keep") %>%
  pivot_wider(
    id_cols = "Method",
    names_from = "Label",
    values_from = "median_LISI"
  ) %>%
  ggplot(aes(x = I(1 - Cluster), y = Batch, color = Method)) +
  geom_point(size = 2) +
  geom_text_repel(aes(label = Method), show.legend = FALSE) +
  labs(x = "1 - cLISI cell type", y = "iLISI batch", )

lisi_plot + labs(title = "LISI Scores")

All Metrics

# arrange the three plots in a single row
prow <- plot_grid(
  kbet_plot + theme(legend.position = "none"),
  asw_plot + theme(legend.position = "none"),
  lisi_plot + theme(legend.position = "none"),
  align = "hv",
  labels = c("A", "B", "C"),
  hjust = -1,
  nrow = 1
)

# extract the legend from one of the plots
legend <- get_legend(
  # create some space to the left of the legend
  kbet_plot + guides(color = guide_legend(nrow = 1)) + theme(legend.position = "bottom")
    
)

# add the legend to the row we made earlier. Give it one-third of 
# the width of one plot (via rel_widths).
plot_grid(prow, legend, ncol = 1, rel_heights = c(1, .1))

Session Information

sessionInfo()

## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] lisi_1.0                    cluster_2.1.4              
##  [3] kBET_0.99.6                 lme4_1.1-31                
##  [5] Matrix_1.5-3                scRNAseq_2.12.0            
##  [7] RSpectra_0.16-1             harmony_0.1.1              
##  [9] Rcpp_1.0.9                  SeuratObject_4.1.3         
## [11] Seurat_4.3.0                r.jive_2.4                 
## [13] gridExtra_2.3               kableExtra_1.3.4           
## [15] scran_1.26.2                scater_1.26.1              
## [17] scuttle_1.8.4               splatter_1.22.1            
## [19] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
## [21] Biobase_2.58.0              GenomicRanges_1.50.2       
## [23] GenomeInfoDb_1.34.9         IRanges_2.32.0             
## [25] S4Vectors_0.36.2            BiocGenerics_0.44.0        
## [27] MatrixGenerics_1.10.0       matrixStats_0.63.0         
## [29] ggrepel_0.9.2               ggthemes_4.2.4             
## [31] ggExtra_0.10.0              cowplot_1.1.1              
## [33] lubridate_1.9.2             forcats_1.0.0              
## [35] stringr_1.5.0               dplyr_1.1.0                
## [37] purrr_1.0.1                 readr_2.1.4                
## [39] tidyr_1.3.0                 tibble_3.1.8               
## [41] ggplot2_3.4.1               tidyverse_2.0.0            
## [43] tictoc_1.1                 
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3                rtracklayer_1.58.0           
##   [3] scattermore_0.8               bit64_4.0.5                  
##   [5] knitr_1.42                    irlba_2.3.5.1                
##   [7] DelayedArray_0.24.0           data.table_1.14.6            
##   [9] AnnotationFilter_1.22.0       KEGGREST_1.38.0              
##  [11] RCurl_1.98-1.10               generics_0.1.3               
##  [13] GenomicFeatures_1.50.4        ScaledMatrix_1.6.0           
##  [15] RSQLite_2.3.0                 RANN_2.6.1                   
##  [17] future_1.31.0                 bit_4.0.5                    
##  [19] tzdb_0.3.0                    spatstat.data_3.0-0          
##  [21] webshot_0.5.4                 xml2_1.3.3                   
##  [23] httpuv_1.6.8                  viridis_0.6.2                
##  [25] xfun_0.37                     hms_1.1.2                    
##  [27] jquerylib_0.1.4               evaluate_0.20                
##  [29] promises_1.2.0.1              restfulr_0.0.15              
##  [31] progress_1.2.2                fansi_1.0.4                  
##  [33] caTools_1.18.2                dbplyr_2.3.1                 
##  [35] igraph_1.3.5                  DBI_1.1.3                    
##  [37] htmlwidgets_1.6.1             spatstat.geom_3.0-5          
##  [39] RcppArmadillo_0.11.4.3.1      ellipsis_0.3.2               
##  [41] backports_1.4.1               bookdown_0.32                
##  [43] biomaRt_2.54.0                deldir_1.0-6                 
##  [45] sparseMatrixStats_1.10.0      vctrs_0.5.2                  
##  [47] ensembldb_2.22.0              ROCR_1.0-11                  
##  [49] abind_1.4-5                   RcppEigen_0.3.3.9.3          
##  [51] cachem_1.0.7                  withr_2.5.0                  
##  [53] progressr_0.13.0              vroom_1.6.1                  
##  [55] checkmate_2.1.0               sctransform_0.3.5            
##  [57] GenomicAlignments_1.34.0      prettyunits_1.1.1            
##  [59] goftest_1.2-3                 svglite_2.1.1                
##  [61] ExperimentHub_2.6.0           lazyeval_0.2.2               
##  [63] crayon_1.5.2                  spatstat.explore_3.0-6       
##  [65] labeling_0.4.2                edgeR_3.40.2                 
##  [67] pkgconfig_2.0.3               ProtGenerics_1.30.0          
##  [69] nlme_3.1-161                  vipor_0.4.5                  
##  [71] rlang_1.0.6                   globals_0.16.2               
##  [73] lifecycle_1.0.3               miniUI_0.1.1.1               
##  [75] filelock_1.0.2                BiocFileCache_2.6.1          
##  [77] rsvd_1.0.5                    AnnotationHub_3.6.0          
##  [79] polyclip_1.10-4               lmtest_0.9-40                
##  [81] boot_1.3-28.1                 zoo_1.8-11                   
##  [83] beeswarm_0.4.0                ggridges_0.5.4               
##  [85] rjson_0.2.21                  png_0.1-8                    
##  [87] viridisLite_0.4.1             bitops_1.0-7                 
##  [89] KernSmooth_2.23-20            Biostrings_2.66.0            
##  [91] blob_1.2.3                    DelayedMatrixStats_1.20.0    
##  [93] parallelly_1.34.0             spatstat.random_3.1-3        
##  [95] beachmat_2.14.0               scales_1.2.1                 
##  [97] memoise_2.0.1                 magrittr_2.0.3               
##  [99] plyr_1.8.8                    ica_1.0-3                    
## [101] gplots_3.1.3                  zlibbioc_1.44.0              
## [103] compiler_4.2.2                BiocIO_1.8.0                 
## [105] dqrng_0.3.0                   RColorBrewer_1.1-3           
## [107] fitdistrplus_1.1-8            Rsamtools_2.14.0             
## [109] cli_3.4.1                     XVector_0.38.0               
## [111] listenv_0.9.0                 patchwork_1.1.2              
## [113] pbapply_1.7-0                 MASS_7.3-58.2                
## [115] tidyselect_1.2.0              stringi_1.7.12               
## [117] highr_0.10                    yaml_2.3.7                   
## [119] BiocSingular_1.14.0           locfit_1.5-9.7               
## [121] sass_0.4.5                    tools_4.2.2                  
## [123] timechange_0.2.0              future.apply_1.10.0          
## [125] parallel_4.2.2                rstudioapi_0.14              
## [127] bluster_1.8.0                 metapod_1.6.0                
## [129] farver_2.1.1                  rmdformats_1.0.4             
## [131] Rtsne_0.16                    digest_0.6.30                
## [133] BiocManager_1.30.20           FNN_1.1.3.1                  
## [135] shiny_1.7.4                   BiocVersion_3.16.0           
## [137] later_1.3.0                   RcppAnnoy_0.0.20             
## [139] httr_1.4.5                    AnnotationDbi_1.60.0         
## [141] colorspace_2.1-0              rvest_1.0.3                  
## [143] XML_3.99-0.13                 tensor_1.5                   
## [145] reticulate_1.27               splines_4.2.2                
## [147] uwot_0.1.14                   statmod_1.5.0                
## [149] spatstat.utils_3.0-1          sp_1.6-0                     
## [151] plotly_4.10.1                 systemfonts_1.0.4            
## [153] xtable_1.8-4                  nloptr_2.0.3                 
## [155] jsonlite_1.8.4                R6_2.5.1                     
## [157] pillar_1.8.1                  htmltools_0.5.4              
## [159] mime_0.12                     minqa_1.2.5                  
## [161] glue_1.6.2                    fastmap_1.1.0                
## [163] BiocParallel_1.32.5           BiocNeighbors_1.16.0         
## [165] interactiveDisplayBase_1.36.0 codetools_0.2-19             
## [167] utf8_1.2.2                    lattice_0.20-45              
## [169] bslib_0.4.2                   spatstat.sparse_3.0-0        
## [171] curl_5.0.0                    ggbeeswarm_0.7.1             
## [173] leiden_0.4.3                  gtools_3.9.4                 
## [175] survival_3.5-0                limma_3.54.2                 
## [177] rmarkdown_2.20                munsell_0.5.0                
## [179] GenomeInfoDbData_1.2.9        reshape2_1.4.4               
## [181] gtable_0.3.1

Dataset 2 Integration

Joey Hastings

2023-03-10

Load Packages

Bacher T-Cell Data

Cell Frequencies by Batch and Cell Cluster

Preprocess Data

Visualize Raw Counts

PVCA

Perform Integration Methods

JIVE

Run Integration Method

Visualization of Integration

kBET

Average Silhouette Width

LISI Scores

Seurat v3

Run Integration Method

Visualization of Integration

kBET

Average Silhouette Width

LISI Scores

Harmony

Run Integration Method

Visualization of Integration

kBET

Average Silhouette Width

LISI Scores

Integration Method Comparison

t-SNE Plots

UMAP Plots

kBET

Average Silhouette Width (ASW)

LISI Scores

All Metrics

Session Information